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Abstract 



We present a quantum kinetic approach for the time-resolved description of many-body effects in photoion- 
ization processes in atoms. The method is based on the non-equilibrium Green functions formalism and 
solves the Keldysh/Kadanoff-Baym equations in second Born approximation. An approximation scheme is 
introduced and discussed, which provides a complete single-particle description of the continuum, while the 
atom is treated fully correlated. 



1. Introduction 

With the development of ultrashort, high-harmonic generated vacuum and extreme ultraviolet (vuv/xuv) 
laser pulses the path towards time- resolved observation of electronic dynamics in plasmas [HEI, a toms and 
condensed matter has been paved, for an overview see e.g. j3]- Recent experiments allow for the investigation 
of electronic motion and relaxation processes on the attosecond (as) time scale [3]. Using as pump-probe 
techniques by combination of an intense femtosecond infrared (ir) and an ultrashort attosecond xuv pulse 
it became possible to directly probe electronic relaxation processes in multi-electron atoms [5] and xuv- 
induced electron shake-up processes in the time domain by means of time-resolved strong field tunneling 
measurements |6j. 

Both experimental scenarios demand for a time-resolved theory of photoionization (PI) in a many-body 
framework. In previous theoretical investigations, powerful tools have been used, such as non-adiabatic 
tunneling theory and single-active electron approaches by means of solving the time-dependent Schrodinger 
equation 0IE]- In a recent paper, these processes were analytically studied in detail (9]. However, all these 
approaches neglect the electron-electron interaction, which may have non-negligible effects, especially in 
the presence of strong laser fields j6]. A promising concept to address these question is multiconfiguration 
Hartree-Fock, e.g. [10]. In this work, we develop an alternative time-dependent many-body approach to 
atomic PI including electronic correlations which is based on non-equilibrium Green functions. 



We aim at describing atomic systems which are initially in equilibrium, i.e. H(t) = Hq for t < to, and 
are disturbed by a time-dependent external potential for t > to- The equilibrium hamiltonian of N electrons 
in the atom is given by (we use atomic units) 
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Figure 1: The Schwinger/Keldysh contour C running from to to t, back to to and to to — i/8 in the complex time plane. "+" (— ) 
marks the (anti-)chronological real-time branch and (|) the imaginary branch. 



where the potential of the nucleus, vfa), as well as the two-particle Coulomb interaction w(ri — r^) = 
|rj — r j | — 1 are assumed to be spin-independent. As we will work in the grand-canonical ensemble, the 
chemical potential contribution [which appears in the density operator] is subtracted for convenience [last 
term in Eq. dIJ]. For times t > to, the atom is disturbed by a time-dependent external field, and the 
hamiltonian is modified: 

JV N N 

H(t) = H + ^2v ext (r i ,t)^^2h(r i ,t) + - 1 ^2w(r i -T j )-fiN. (2) 

i—l i—1 i^j 

With the last equation we defined the total single particle hamiltonian, h = ho + v ext , which will be used 
below. In this paper we consider the perturbation by an electromagnetic wave in dipole approximation, 
i.e. the field is assumed homogeneous on the scale of the atom, v ext (r,i) = — eE(f) • r. The electric field 
envelope E(i) is assumed to have a Gaussian shape 

£(t) = £ cos[uj(t - t m id)] exp 

with a pulse duration r. 



(t — tm id) 

2r 2 



<2 1 



(3) 



2.1. Contour Green functions 

In the Keldysh/Kadanoff-Baym approach, the central quantity is the one-particle nonequilibrium Green 
function G, which is the time-ordered expectation value of the product of two field operators: 

G(l,2) = -i(T c [*tf(l)fck(2)]>, (4) 

where the variable 1 = (ri,<7i,ii) comprises position, spin projection and time, and the field operators are 
considered in Heisenberg representation. In the following we denote X! = (ri,<Ti). The ensemble average in 
([4} is performed in the grand-canonical ensemble, i.e. with the trace over the unperturbed grand canonical 
density operator Z~~ 1 e~ l3H ° . The Green function is defined on the Schwinger/Keldysh time contour C, see 
Fig. |2.1| which allows for an extension of the groundstate and equilibrium formalism and the diagram 
technique, to non-equilibrium jTT], for an overview see |12l 113] . The contour runs from the initial time to to 
the current time t, back to to and, finally to to—i(3 in the complex plane where j3 — 1/fcflT. The time-ordering 
operator Tq in equation Q arranges operators with time arguments "later" on the contour to the left. The 
propagation along the complex branch corresponds to the Matsubara formalism, in which the equilibrium 
density operator is expressed by a time-evolution operator in complex time, e~P Ha = U(to — ij3, t ). With 
this, the Green function ([3} takes the explicit form 

G(1 2) = Tr{U(to-ip,t )T e [* H (l)¥ H (2)}} 
Tr{U(t - i[3,t )} 
2 



in which the time-arguments t\ and t% each lie on one of the three branches of the contour C. It is now 
convenient to introduce a set of subordinated Green functions, which depend on the location of the two 
time-arguments on the contour. Accordingly, the Green function becomes a 3 x 3 matrix, 

G c G< \ /G++ G+- G+l\ 

G> G a := G~+ G G~\ , (6) 

G^ G^ G M j \G\+ G\~ Gil/ 

where {+, — , |} mark the position of the respective time-argument on the branches, where the left (right) 
symbol corresponds to the first (second) time argument. For the notation in the right part of Eq. see 
Fig. l.The notations in the left part of Eq. ^ show the relation to the standard definitions, where G^ 
denote the correlation functions with two real-time arguments and G^ [G^] denotes the correlation function 
in which the first [second] argument lies on the imaginary branch and the second [first] on one of the real 
branches. Finally, G M denotes the Matsubara (imaginary time) Green function of equilibrium theory. The 
matrix notation allows one to eliminate the time contour C and to consider, in the following, only functions 
of real time arguments. 

Of the four real-time functions {G c ,G < ,G a ,G > }, only two are linearly independent, which is why we 
consider in the following only the correlation functions G < and G > . The mixed Green functions G^ and 
G^ account for the evolution of the initial equilibrium state, which itself is determined by the Matsubara 
Green function G M . As the equilibrium Green function only depends on the difference of two complex 
time-arguments, we consider the real function G M (x 1 ,x 2 ,r 1 — t 2 ) := — iG(xi — £ti,x 2 — ir 2 ), which is 
defined in the interval [— /3, (3\ and obeys the symmetry G M (xi, x 2 , t) = — G M (xi, x 2 , t — j3). It is, therefore, 
sufficient to determine G M in the range [— /3, 0]. 

2.2. Keldysh/Kadanoff-Baym equations 

To compute the time evolution of a multielectron atom, we need to solve the equations of motion of 
the Green function G. The equations for the Keldysh matrix function are the Keldysh/Kadanoff-Baym 
equations (KKBE) defined on the contour C |14] . 

{id tl -h(l)}G(l,2) = S c (l-2)+ f d3£[G](l,3)G(3,2), (7) 

Jc 

which is to be supplemented by the corresponding adjoint equation. S[G] denotes the irreducible self-energy 
which is also a 3 x 3 matrix containing mean field (Hartree-Fock) and correlation effects. We will discuss 
approximations to the self-energy in section |2.3[ 

The matrix equation Q is equivalent to a coupled system of equations for the subordinated Green 
functions which are derived applying Langreth's rules to the right hand side of Eq. Q, e.g. [12], and have 



the form (for compactness we suppress space and spin variables) 

{d Tl - h a } G M ( n - r 2 ) = 5(n - r 2 ) + f df S M (n - f) G M (f - r 2 ) , (8) 

Jo 

{id tl -h(h)} G > (t 1 ,t 2 ) = 7 > (ti,t 2 ), (9) 

{-id t2 -h(t 2 )} G < (t 1 ,t 2 ) = ^{tuh), (10) 

{id tl -h(h)} G\h,T 2 ) - /l(ti,T 2 ), (11) 

{~id t2 -h(t 2 )} G r (n,< 2 ) = J r (Ti,i 2 ), (12) 



Complex time arguments are indicated by r and the system has to be supplemented by the adjoint equations. 
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Figure 2: Used approximation for the self-energy, from left to right: Hartree and Fock term, second Born and second order 
exchange contributions. 

The collision integrals are given by (without loss of generality, from now on we use to = 0) 



dt 



o 



Z R (t 1 ,t)G>(i,t 2 ) + i:>(t 1 ,t)G A (t,t 2 )\ -i df El(ti,r)Gr(r,t 2 ), (13) 



Jo 



G fl (t 1 ,t)E<(f,t 2 ) + G > (ti,t)S A (t,t 2 ) -i/ dfGHti.fjS^f.ta), (14) 

J Jo 

■P 

(15) 



/ 1 (tl,T 2 )= f 2 dt^ R (t 1 J)G\t,T 2 )+ f df^(t 1 ,f)G M (f-T 2 ), 
J t JO 



and I^(ti ,i 2 ) = /I (ti, r 2 +/?). The last terms in Eqs. ( 13"fl5 l account for the evolution of initial correlations. 
To shorten the notation, here we have introduced the retarded (R) and advanced (A) Green functions which 
are defined according to G R ^ A {t\,t 2 ) — ±0[±(ii — t 2 )][G > (t\, i 2 ) — G < (ti,t 2 )] and analogously for the 
self-energies. 

The equations of motion have to be supplemented by initial conditions. In this paper we start from 
the correlated equilibrium state of the electrons in the atom which is determined by the Matsubara Green 
function G M - the solution of the equilibrium Dyson equation, Eq. (|8}. For a convenient numerical treatment, 
it is transformed into the integral form [15] (we use r = t\ — r 2 ): 



G m (t) = G°(r)+ fdfdf G°(r-f) E[G M ](f-f) G M (f). 



(16) 



G° is a reference Green function - in our case given by the Hartree-Fock Green function - and E(r) = 
E m (t) — (5(t)£°, where E° is the reference self-energy evaluated with G°. The boundary conditions for the 
Dyson equation (16 1 are given by the Kubo-Martin-Schwinger condition, G M (r) = G M (r — i/3). 



2.3. Self-energy approximations 

The equations of motion for the Green functions are formally exact if the self-energy would be known. 
Here we use a conserving approximation [16] for E which includes Hartree-Fock contributions and correlations 
in the second Born approximation with exchange. The corresponding Feynman diagrams are shown in Fig. 
2.3 and correspond to the explicit expressions (defined on the time contour) 



E HF (1,2) = iG(l,2)to(l,2) - i<y c (l-2) J d3to(l,3) G<(3,3) , 

E 2B (l,2) = i 2 f d3d4G(l,3)tu(l,4)G(3,4)G(4,2)ti;(3,2) - i 2 [ d3d4G(l,2)w(l,3)w(2,4)G(4,3)G(3,4) , 
Jc Jc 

where we denoted w(l,2) = 6c(l — 2)w(ri — r 2 ). Obviously, the Hartree-Fock self-energy is time local. 

The advantage of using the method of nonequilibrium Green functions is that it provides a fully selfcon- 
sistent approach to electronic correlations in atoms in equilibrium and nonequilibrium. Solving the KKBE 
with the above self-energies one obtains the time evolution of a many-body system, thereby fully preserving 
momentum, angular momentum and total energy. Furthermore, due to the inclusion of memory effects 
(time integrations in the collision terms) no restriction with respect to the times apply, which is particularly 

4 



important for ultrafast processes in optically excited atoms. Finally, the electromagnetic field is included 
non-perturbatively which allows to investigate the nonlinear dynamics of atoms in the presence of a strong 
excitation. 



3. Implementation 

3.1. Basis representation 

Despite their mentioned above attractive properties, the KKBE 
Already for one-dimensional systems, they constitute a set of four-dimensional integro-differential equations 
(not counting spin degrees), whereas for three-dimensional systems they are even eight-dimensional. Without 
further approximations, this is far beyond today's numerical possibilities. A first way around this is to 
expand all quantities in terms of suitable single-particle basis functions {4>k} with k = 1,2,.. . JVb as was 
demonstrated in Ref. [T7]. For example, the Green function and, likewise any other single-particle quantity, 
becomes an Nb x Nb dimensional matrix with the elements 

Gyfa.ia) = ydxidx 2 0*(xi)G(l,2)^(x 2 ), i,j = l,...N b . (17) 

Below we will use for {4>k} a set of Nb orthonormal Hartree-Fock orbitals constructed from atomic orbitals. 
The expansion of the two-particle interaction term yields a tensor with four indices, the two-electron integrals: 

Wijkl = J dxidn 2 4>l(xi)4>j(xi) w(xi,x 2 ) </>fc(x 2 )0i(x 2 ) , (18) 

which are time-independent. The main advantage of the expansion is that we have eliminated the coordinate 
dependencies. In terms of the basis, all equations considered earlier become equations for matrices depending 
on two time arguments. 

So far, we have not explicitly specified the single-particle basis, the results are completely general. 
Also, the treatment of systems with different dimension is conceptually greatly simplified since it requires 
nothing more than a respective set of electron integrals, i.e. the matrix representation of the single-particle 
hamiltonian h and its parts h pot ,hjji n , the overlap matrix Oij = J dx.(f>*(x)<pj(x) as well as the electron 
repulsion integrals Wijki- Furthermore, for the laser excitation the dipole matrix d = — er is needed. Up 
to now, we have implemented one-dimensional numerical orbitals, three-dimensional numerical orbitals for 
central potentials, Slater type orbitals for atoms and Gaussian type orbitals for arbitrary molecules. Within 
this work we will illustrate the method for a one-dimensional model atom. 



9JT2 I are very hard to solve numerically. 



3.2. Solution procedure 

We briefly outline of the solution procedure. We follow the techniques developed earlier, see e.g. |18l [T9l 
IT?] and references therein, a detailed description of the algorithm is given in Refs. |15] and |20] . 



1. A single-particle basis {<j>k} is chosen, and the one- and two-electron integrals (18 1 are calculated. 

2. This provides the input for a Roothaan-Hartree-Fock calculation [21] yielding the HF energies and 
orbitals. 

3. The electron integrals are transformed to the HF basis [22], and the reference Green function is set up 
on a uniform power mesh, an adapted imaginary time-grid, as G^-(t) = Sij ni e~( Ci ~^) r , for t £ [— /3, 0). 
rii is the occupation number of the HF-orbital i, which is determined by the Fermi distribution. 

4. The Dyson equation ([8} is solved iteratively, until a self-consistent solution for the correlated Matsubara 
Green function G M is found. Thus the equilibrium problem is solved. 

5. The function G M determines the initial conditions of the time propagation given by the KKBE, Eqs. 
(9][l2l according to 

G>(0,0) = tG M (Q-), G>(0,0) = tG M (0+) = -iG M (-/3), (19) 

Gl(0,-«r) = iG m {-t), G r (-zr,0) = iG m {t) = -iG m {t-[5). (20) 
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Figure 3: (a) Schematic view of the ID model Beryllium atom. Shown is the square of the lowest two doubly occupied 
Hartree-Fock orbitals (HFO) and the ground state density, (b) Density of states (DOS) for N b = 100 Hartree-Fock basis 
functions obtained after convolution with a Gaussian of FWHM = 0.02. Grey curve is the result for 100 HF basis functions 
and additional 200 ideal basis function to better resolve the continuum, see Sec. [5] 



6. The KKBE are rewritten in terms of the time-evolution operator U(t+At, t) = Q- lh ( t ) At and are solved 
in the two-time plane by standard techniques for ordinary differential equations. We have currently 
implemented a fourth-order Runge-Kutta scheme [23]. Due to the symmetry G<(ti,t2) = [G<{t2,ti)\ * 
and the boundary condition G > (t,t) = G < (t,t) — i, it is sufficient to propagate the lesser Green 
function in the upper triangle t\ > £2 of the real-time plane and the greater Green function in the 
lower triangle t\ < £2- The mixed Green functions are both propagated in the whole real-complex 
time-plane, though, in the long-time limit the initial correlations decay. 



4. Simulation results 



4-1. One-dimensional model atom 

To examine the presented formalism, we consider a one-dimensional model atom defined by a regularized 
Coulomb potential, v(x) — —Z(x 2 + k^) -1 / 2 , where Z is the atomic number. The screening-parameter 
k_ is introduced to avoid the computational difficulties arising from the singularity at the origin. Likewise, 
the two-particle interaction is modified according to w(xt,X2) = [{xi — X2) 2 + ft+] -1 / 2 • This model has 
been used successfully in many studies of atom-laser interaction where also the influence of the choice of the 
screening parameters has been investigated. Here we follow Ref. |2D and use k_ = k + = 1. 

We illustrate the method by considering beryllium (N = Z = 4). A sketch of the confinement and the 
HF basis functions shifted by the orbital energies is shown in Fig. [4] (a). As can be seen, there are two 
bound states at the energies E$ — —1.371 and E\ = —0.312, each occupied by two electrons. In Fig. [4] (b) 
the density of states (DOS) is plotted, which is obtained by a modeling the continuum by a box of width 
200 a.u. and after convolution of the delta-peaks with a Gaussian of FWHM = 0.02. It shows the typical 
l/V^-decay. 



(.2. Ground state properties 

Results for the ground state energies computed in HF and second Born approximation are presented 



in Table 4.1 and compared to an exact diagonalization (CI) calculation. For the present comparison it is 
sufficient to consider a CI calculation using a basis of size Nb = 15 which is performed sufficiently fast. For 
the solution of the Dyson equation, a maximum number of Nb = 80 basis functions has been used. Together 
with the scaling properties, which are barely affected by the particle number but mostly by the number of 
basis functions, we are able to go far beyond the region accessible by sophisticated (full) CI methods |25] . 
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HF 


2ndBorn 


CI 


15 


-6.7390 


-6.7694 


-6.7831 


30 


-6.7393 


-6.7706 




80 


-6.7395 


-6.7710 





Table 1: Ground state energies of the ID Be model atom obtained from Hartree-Fock, second Born and Configuration Interaction 
calculations using different basis dimensions JV&. 



To obtain the ground state within the finite temperature formalism, an inverse temperature of /3 — 100 is 
used. 

As can be seen, already the HF energies are close to the CI results. The inclusion of correlations on the 
second Born level yields a further improvement, accounting for 69 % of the correlation energy. This confirms 
the trend observed for the ground states of real atoms [20 . In this reference, it is also demonstrated 
how to obtain the ionization energies from the Green function using the Extended Koopmans theorem. 
A corresponding calculation for the one-dimensional Beryllium model yields a first ionization energy of 
I p = 0.303, which should be compared to the Hartree-Fock ionization potential from the (conventional) 
Koopmans theorem, I p = 0.312, which is known to overestimate the ionization energy. 

4-3. Time- dependent ionization dynamics following a short UV pulse 

Let us now consider the perturbation of the atom by different electromagnetic pulses. We use a Gaussian 
pulse (|3| with a fixed number of cycles (r = 107r/w) and an amplitude £ = 0.1. We consider three 
frequencies: uj\ = |£aJ/2, L02 = {\Eq\ + \Ei\)/2 and W3 = 1.2 • \Eq\. A classification of the different laser 
pulses is given by the Keldysh-parameter |26] 7 = <J I p /2U P , where I p is the ionization potential, and the 
ponderomotive potential is given by U p = £q/Auj 2 . The system is propagated to T — 200 in Hartree- 
Fock approximation with Nb = 100 HF basis functions. Fig. |4.3| shows the time-dependent results for 
the level populations during the pulse, as well as the time-dependent expectation values of the electron 
coordinate (ar). In the upper pictures (oj\ — \Ei\/2 = 0.156) the frequency is too small to allow for PI of an 
electron. Nevertheless, strong ionization of the upper electrons is observed which is due to tunnel ionization 



(71 = 1.25). The expectation value of the electron position, Fig. 4.3 (b), confirms this interpretation. For 
the frequency u>2 = 0.842 PI of the upper electrons and, for the case cu^ — 1.645, direct PI of electrons from 
both levels is possible. In these cases, the Keldysh parameters belonging to the second orbital is 72 = 6.7 
and 73 = 13, respectively, thus tunneling is not relevant, but multi-photon ionization occurs. 

We have also performed various correlated simulations of the PI dynamics using the second Born approx- 
imation for the selfenergy. However, they require a propagation of the Greens functions in the full two-time 



plane by solving the KKBE, making the calculation computationally very costly. The upper row of Fig. 4.3 
also contains correlated results for the dynamics up to t = 25. During the initial phase the deviations form 
TDHF are still small. 



5. Approximate treatment of continuum states 

The main restriction in the time-dependent description of PI on the level of the Born approximation 



is the very large CPU time and memory requirement, the result shown in Fig. 4.3 was obtained using 300 
time steps with Nb = 20 basis functions. The calculations took 48 hours on a single CPU and required 
4 GB of main memory. For the first part - the equilibrium calculations - the basis representation saves 
a lot of numerical effort compared to a solution in coordinate representation. HF orbitals turn out to 
provide an efficient basis also for the correlated equilibrium Green function, allowing to restrict the basis to 
a dimension of the order Nb < 100 since normally only the lowest orbitals give the dominant occupation. In 
nonequilibrium, however, it is generally not possible to truncate the basis since, in principal, every orbital 
can be occupied during the excitation process. In particular, PI will lead to occupation of continuum states, 
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Time (atomic units) Time (atomic units) 



Figure 4: Dynamics of the system induced by three different pulses with electric field strength E = 0.1 and three different 
frequencies in Hartree-Fock (HF) and second Born (2B Approximation. For HF a basis of dimension N b = 100 was used. 
For 2B N b = 20 and the system was propagated up to T = 25 a.u. (a) Time-dependent occupation numbers of the occupied 
Hartree-Fock orbitals, the upper (lower) curve refers to the first (second) orbital, (b) Time-dependent expectation value of 
position. The insets in the upper plots show the comparison between HF and 2B more in detail. 



and at high intensities, the electron energies may become very large [27J. Obviously, a fully correlated 
description of the whole continuum is prohibitive. At the same time, an electron "born" in the continuum 
with large kinetic energy will be only very slightly disturbed by correlations with the other electrons. This 
naturally suggests to develop an approximation scheme which is based on a sub-division of the basis into 
low and high lying orbitals which are treated with different levels of accuracy with respect to many-body 
effects. As a result we may hope that the ionization dynamics of the atom can be resolved sufficiently well, 
both in energy and in time, including electronic correlations with the required accuracy. 

Our ansatz is the following: We divide the basis {</>&} in three sub-systems, Nb = N C011 ' + N HF + N ld . 
The first sub-system (1) contains the N corr energetically lowest orbitals (for example the lowest atomic 
bound states) and is treated fully correlated. The second (2) contains a number of N HF low-lying continuum 
orbitals which are treated in Hartree-Fock approximation whereas the third (3) of dimension N ld is treated 
without any particle-particle interaction. In all three sub-systems the single-particle contributions, such as 
the external electromagnetic field are included exactly, thus fully taking into account non-linear effects. For 
a given Nb, the sub-division of the basis and the numbers N colr , N HF of functions in the basis parts are 
arbitrary and can be adapted to the considered atom and excitation conditions. 

Let us now introduce this sub-division into the nonequilibrium Green functions scheme. All Green 
functions (all Keldysh components) and the single particle hamiltonian now become 3x3 matrices, 

(Gn G12 Gi3\ /hn hi2 hi3\ 

Gai G 22 G 23 , h = h 21 h 22 h 23 . (21) 
G31 G32 G33/ \li31 h.3 2 1133/ 



The self-energy which, in the following, is separated into a time-diagonal Hartree-Fock and a two-time 
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correlation part, is intrinsically of the same structure. Following the idea of our approach, approximations 
will be introduced by systematically neglecting certain blocks of the self-energy matrix: 



corr 



























" m - I £«! £° 2 ] , S corr = ( 0], (22) 

0y 

where S° denotes the HF selfenergy evaluated with uncorrelated Greens functions. With this, the KKBE 
attain the form 



{ia 4l -h(t 1 )}G(t 1 ,t 2 ) =<fc(l-2) 





^ S corr Gii 




£COrr 


dt 3 













V o 









G 
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(*i,*a), (23) 



where the HF self-energy has been included in the mean-field hamiltonian h(ti) := h(ti) + £ HF (fi) on the 
l.h.s. With these notations, Gn contains the Green functions of electrons occupying atomic bound states or 
undergoing transitions between low lying bound states whereas the information about ionization processes 
is contained in G12 and G13. Electrons in the continuum are described by G22 and G33. This model is 
closely related to the Bloch equations of atomic physics or their generalizations to semiconductor optics, for 
a formulation using nonequilibrium Green functions, see e.g. |12l 128] . A similar scheme has recently been 
reported in Refs. |29 j and [30 , where it was applied to quantum transport. 

If the ionization is weak, as is normally the case with an (X)UV pulse produced from an optical laser 
via high harmonics, the ionization components G12 and G13 will be much smaller than contributions from 
occupied orbitals in Gu. We then may expect that correlation effects in Gi 2 and G13 play a minor role. 
This allows us to further simplify the model by setting S corr Gi2 ~ S corr Gi3 ~ on the right hand 



side of Eq. (23). This, obviously, becomes questionable in the case of intense fields, such as optical or IR 
probe beams etc., but this question is beyond the scope of this work and will be considered in a forthcoming 
analysis. With this approximation, only the Green function of system (1) has to be considered correlated 
and only Gu has to be evolved on the full two-time plane whereas the other Green functions are completely 
determined by the information on the time-diagonal. 

With these approximations we can write down the final system of equations to be solved as 

3 - r 
i^Gufa.taJ-VhyfrjG^ti.ta) = S c (t l7 t 2 ) + / dt 3 S corr (^, t 3 )G n (t 3 , t 2 ) , 

Jc 

3 

id tl G tk {t 1 M)~Y,^( t i) G ok{hM) = 6c(ti,t 2 ), (*,*) = {(1,2), (2,1), (2,2)}, 
3 



idt 1 G ik (t 1 ,t 2 )-J2h ij (t 1 )G jk (t 1 ,t 2 ) = Sc(t 1 ,t 2 ), i or k = 3. (24) 

We give some further information on the performance of the approximation. The limiting factor in the 
formalism is the number of basis functions. Based on our first tests [3T we expect that our scheme is 
capable to treat correlated subsystems of a size of the order of N corr = 50. This is sufficient to describe 
the nonequilibrium ionization dynamics of many atoms and small molecules. The size of the Hartree-Fock 
subsystem can be substantially larger since it depends only on a single time-argument. There the limiting 
factor is the computational effort of the two-electron integrals, which grows as 0[(jV hf ) 4 ]. Finally, the size of 
the ideal sub-system can, in principle, be chosen nearly as large as in standard solutions of the single-particle 
time-dependent Schrodinger equation. 

Here, we present the first equilibrium results of our approximation scheme. One example is seen in 
Fig. I^Jd. where, in addition to the DOS calculated for a pure HF basis with N HF = 100, we include the 
results for a basis supplemented by N ld = 200 ideal basis functions (N corr = 0). The figure shows that 
this allows us to extend the basis to significantly higher energies and better resolve the continuum. This 
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Figure 5: (a) Ground state energies obtained with the approximation scheme (open symbols), full HF (crosses), full CI (filled 
squares) and full second Born calculations (filled triangles) versus basis size. For the full calculations the abscissa indicates the 
number of basis functions. For the approximation scheme the abscissa denotes the number N corr of correlated orbitals, which 
are extracted from a Hartree-Fock calculation with Nb = 15. (b) DOS for a HF calculation with 10 basis functions, (c) DOS 
for a HF calculation with 15 basis function (light grey) and the approximation scheme using 

N corr _ 1Q an( j N HF _ g Dark 

area show the DOS corresponding to the correlated orbitals. 



opens the way to investigate significantly higher excitations than before, which is particularly important 
for PI with (x)uv photons or for multiphoton ionization processes. The second example is shown in Fig. 
[5] where we compute the ground state energy of the model atom for varying basis subdivisions, N corr and 
N HF at fixed total basis size N corr + N HF = N b = 15 (for simplicity, iV id = 0). The results are compared 
with independent solutions of the Dyson equation and CI calculations for different basis sizes Nb without 
applying the approximation scheme. While the full calculations show fast convergence, with increasing Nj,, 
the approximation scheme practically remains at the HF level for N corl < 10. Only for N corl > 10 the 
results become better than the HF ones and eventually approach the correlated values. This behavior is not 
a property of the second Born approximation but is also observed in the CI calculations. 

This slow convergence is unexpected, since one would anticipate a larger HF basis to be more adequate in 
describing correlation effects. To analyze the reasons of this behavior we show in Fig.[5}3) and c) the DOS for 
two HF basis dimensions N HF = 10 and N HF = 15, respectively. As one can see, the larger basis causes an 
upshift of the two highest peaks of the DOS and a strong increase of the central peak around E = 0. While 
this has only little effect on the HF ground state energy it strongly influences the convergence behavior 
of our approximation scheme. For illustration, in Fig. |5]c) we also show the DOS for a calculation with 
N corr = 10 and N HF = 5 (dark area) which is close to the DOS for a pure HF calculation with N HF = 15. 
Only when N colr is so large that the correlated orbitals extend beyond the high peak is the approximation 
scheme approaching the correlated calculation. We are presently investigating how to avoid this unwanted 
behavior of our scheme in order to achieve a faster convergence. This will also be the basis for extending 
this approximation scheme to nonequilibrium calculations of PI processes. 



6. Conclusion and outlook 

In this paper we applied the nonequilibrium Green functions approach to atomic photoionization. Due 
to the large basis size needed for the description of the continuum, which make fully correlated calculations 
unfeasible, we have derived an approximation scheme. Starting from the exact Keldysh-Kadanoff/Baym 
equations for the correlation functions we neglected certain matrix elements in the basis representation 
and obtained a scheme which is expected to be efficient for application to atoms. The idea is to restrict 
the computationally costly evaluation of the correlation self-energy to a subset of the single-particle basis. 
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Thereby only the lowest lying bound states are treated fully correlated, while the continuum is approximated 
by Hartree-Fock and ideal basis functions. We have performed the first tests on the ID-Beryllium model 
atom which provide a framework for further investigations. Further optimization of the basis subdivision 
towards improved convergence and application to nonequilibrium photoionizaiton dynamics will be presented 
in a forthcoming work. 
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